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Abstract 

The theory of non-oscillatory scalar schemes is developed in this paper in terms of the local extremum dimin- 
ishing (LED) principle that maxima should not increase and minima should not decrease. This principle can 
be used for multi-dimensional problems on both structured and unstructured meshes, while it is equivalent to 
the total variation diminishing (TVD) principle for one-dimensional problems. A new formulation of symmet- 
ric limited positive (SLIP) schemes is presented, which can be generalized to produce schemes with arbitrary 
high order of accuracy in regions where the solution contains no extrema, and which can also be implemented 
on multi-dimensional unstructured meshes. Systems of equations lead to waves traveling with distinct speeds 
and possibly in opposite directions. Alternative treatments using characteristic splitting and scalar diffusive 
fluxes are examined, together with a modification of the scalar diffusion through the addition of pressure 
differences to the momentum equations to produce full upwinding in supersonic flow. This convective upwind 
and split pressure (CUSP) scheme exhibits very rapid convergence in multigrid calculations of transonic flow, 
and provides excellent shock resolution at very high Mach numbers. 


1 Introduction 

Over the past decade the principles underlying the design of non-oscillatory discretization schemes for com- 
pressible flows have been quite well established. A very large number of variations of artificial diffusion, 
upwind biasing and flux splitting have been proposed and tested [22, 38, 26, 34, 31, 11, 40, 17]. In the same 
period multigrid acceleration schemes have also been the subject of widespread investigation, and have proved 
effective, particularly for subsonic and transonic flow. The use of limiters to enforce monotonic solutions, has 
proved, however, to have an adverse effect on multigrid convergence. If fact, it is not uncommon for schemes 
with limiters to become trapped in limit cycles. Limiters also tend to reduce the accuracy of solutions, par- 
ticularly in regions containing smooth extrema. The first purpose of this paper is to develop a systematic 
procedure for the analysis and design of a broad class of schemes which satisfy monotonicity constraints on 
both structured and unstructured grids, and to illuminate some of the connections between alternative formu- 
lations. The second purpose is to show ways in which these schemes can be modified to improve both their 
accuracy and their rate of convergence to a steady state. Schemes which blend low and high order diffusion 
[22], and both symmetric and upstream constructions using anti-diffusive terms controlled by limiters [19], 
are readily included within the framework of this paper. The connection between schemes of this type and 
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schemes which require the exact solution of a Riemann problem at each cell interface has been widely examined 
elsewhere [12, 9, 45], and is not explored here. 

Two main issues arise in the design of non-oscillatory discrete schemes. First there is the issue of how to 
construct an approximation to a scalar convection or convection- diffusion equation which is non-oscillatory, 
captures discontinuities with high resolution, and is sufficiently accurate. Second there is the issue of how to 
construct a numerical flux for a system of equations with waves traveling at different speeds, and sometimes in 
opposite directions. These two issues can be treated essentially independently, and by combining alternative 
non-oscillatory formulations with different constructions of the numerical flux one arrives at a matrix of 
candidate high resolution schemes, all of which may have acceptable characteristics. Reference [42] examines 
the performance of such a matrix of schemes for viscous boundary layers. 

Section 2 reviews the conditions for the construction of non-oscillatory schemes for scalar conservation laws. 
Following a line adhered to in a number of works [6, 47, 23, 37] , including several by the present author 
[16, 17, 20], it is suggested that the principle of non-increasing maxima and non-decreasing minima provides 
a convenient criterion for the design of non-oscillatory schemes. This principle contains the concept of total 
variation diminishing (TVD) schemes for one-dimensional problems, but can readily be applied to multi- 
dimensional problems with both structured and unstructured grids. Such local extremum diminishing (LED) 
schemes can be realized by making sure that the coefficients of the discrete approximation are non-negative. 
First order accurate schemes satisfying this principle are easily constructed, but are too diffusive. It is well 
known that schemes which strictly satisfy the LED principle fall back to first order accuracy at extrema even 
when they realize higher order accuracy elsewhere. This difficulty can be circumvented by relaxing the LED 
requirement. Therefore the concept of essentially local extremum diminishing (ELED) schemes is introduced. 
These are schemes for which, in the limit as the mesh width Ax — ► 0, maxima are non-increasing and minima 
are non-decreasing. 

One approach to the construction of high resolution schemes which combine monotonicity and higher order 
accuracy is to blend low and high order diffusive terms as, for example, in the Jameson-Schmidt-Turkel (JST) 
scheme [22]. It is proved in section 2.2 that with appropriately chosen coefficients the JST scheme is LED. 
Moreover, the coefficients can be chosen so that the scheme is both second order accurate at smooth extrema 
and ELED. Another approach to the construction of higher order schemes, which has been adopted by several 
authors [47, 15, 44, 46] is to add limited anti-diffusive terms to a lower order scheme. In section 2.3 and 
2.4 this procedure is used to derive a general family of symmetric limited positive (SLIP) schemes for both 
structured and unstructured meshes. Moreover, the switch between low and high order terms in the JST 
scheme can be formulated in such a way that the JST scheme becomes a special case of the SLIP scheme. 
A slight modification of the SLIP formulation produces a corresponding family of upstream limited positive 
(USLIP) schemes, which are derived in section 2.5, and resemble some well known upwind schemes [30, 40]. 
The limiters in the SLIP and USLIP schemes can be relaxed so that the schemes are second order accurate 
at smooth extrema and ELED. Section 2.6 shows how a sequence of successively higher order ELED schemes 
can be derived. 

Section 3 discusses the treatment of systems of equations with several dependent variables. In order to apply 
the local extremum diminishing (LED) principle, the flux may be split in a manner which corresponds to the 
characteristic fields, so that the scheme is designed to limit extrema of the characteristic variables. The Roe 
flux [34] provides a way to produce schemes that resolve stationary shock waves with a single interior point. 
The use of a scalar diffusive flux constructed directly from the solution variables leads to simpler schemes 
which can resolve shock waves with several interior points, and exhibit no overshoots provided that enough 
diffusion is introduced locally. These schemes have proved quite effective for steady state calculations. Very 
rapid convergence to a steady state can be achieved by the introduction of multigrid acceleration techniques. 
Because of their low computational costs scalar diffusive schemes have proved quite suitable for industrial use, 
and they have been successfully used for aerodynamic analysis in the design of aircraft such as the YF-23 [7]. 

Scalar diffusion has the drawback that in order to stabilize the calculation, it tends to introduce more 
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diffusion than is really needed. In the treatment of high Reynolds number viscous flows in which the viscous 
effects are mainly confined to thin boundary layers, it is important to keep the numerical diffusion in the 
boundary layers as small as possible. In supersonic flow the region of dependence is purely upstream, while 
the use of scalar diffusion cannot produce a discrete scheme which is fully upwind. This can be remedied by 
the introduction of pressure differences in the momentum equation. It is then possible to construct a first order 
scheme which essentially reduces to pure upwinding in supersonic flow, and which may be used as the basis for 
constructing higher order schemes. Upwinding of the pressure requires the introduction of terms which may 
be related to the wave particle scheme of Deshpande, Rao and Balakrishnan [33, 4], and flux splitting recently 
proposed by Liou and Steffen [27], This simple convective upwind and split pressure (CUSP) scheme produces 
discrete normal shock waves which contain two or three interior points in transonic flow, and become sharper 
at very high Mach numbers. 

Section 4 presents results of various test calculations for one-, two- and three-dimensional problems. It is 
verified that the JST and SLIP schemes with Roe flux provides high resolution of shock waves in shock tube 
simulations. They also produce shock waves with one interior point in steady transonic flow calculations for 
airfoils. The CUSP scheme has the advantage, however not only of reduced computational complexity, but also 
of significantly faster convergence to a steady state. Transonic solutions using the CUSP scheme on a 160x32 
grid are presented for three different airfoils. Each was obtained in 12 multigrid W-cycles with a multistage 
explicit time stepping scheme. The reduction of the number of steps needed for global convergence to 12 is 
the culmination of 12 years of effort. 


2 Non-oscillatory schemes for scalar equations 

2.1 Local extremum diminishing (LED) and essentially local extremum dimin- 
ishing (ELED) schemes 

Consider the discretization of a time dependent conservation law such as 

dv d x d , . /1X 

a + to /< ' ,) + V w = 0 ' () 


for a scalar dependent variable v 
points are numbered in some way, 
is expressed in semi-discrete form 


Then on introducing Taylor series 
term 


on an arbitrary (possibly unstructured) mesh. Assuming that the mesh 
let Vj be the value at mesh point j. Suppose that the approximation to (1) 
as 


dvj 
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expansions for v (x* — Xj , t/* — Vj) 


it follows that in the absence of a source 


E c j* = 0 

k 

Thus there is no loss of generality in writing the scheme as 


dvj 

Hi 


= '£,c jk (v k -v i ). 
k*i 


Suppose that the coefficients are non-negative 


Cjk >0, k ^ j. (2) 

Then the scheme is stable in the norm, since if vj is a maximum, vt-Vj < 0, so that ^ < 0, and similarly 
a minimum cannot decrease. Suppose, moreover, that the stencil of the discrete scheme is compact 


Cjk = 0 if j and k are not nearest neighbors 


(3) 



Then if Vj is a local maximum (over the stencil of the difference scheme) Vk ~ Vj <0, with the consequence 


that ^ < 0. Thus a local maximum cannot increase, and similarly a local minimum cannot decrease. Such 
a scheme will be called local extremum diminishing (LED). 

This criterion has been proposed by various authors [6, 17, 21, 23, 47] as a convenient basis for the construc- 
tion of non-oscillatory schemes on both structured and unstructured meshes. It assures positivity, because if 
v is everywhere positive, then its global minimum is positive, and this cannot decrease. When specialized to 
one dimension it also leads to the class of total variation diminishing (TVD) schemes proposed by Harten [11]. 
The total variation of v is 


tv( v) = r 

J — 0 


dv 

dx 


dx } 


that is the sum of the absolute values of the variation over each upward and downward segment. It was 
observed by Laney and Caughey [23] that each extremum appears in the variation of the segment on each side 
of that extremum, with the consequence that 


TV (v) — 2 maxima — mmirnaj , 


if the end values are fixed. Thus, if a one-dimensional scheme is LED, it is also TVD. On a triangular mesh, 
a definition of total variation such as 


TV (u) = J ||Vti|| dS 


is not an entirely satisfactory measure of oscillation. This is illustrated in Figure 1, where the total variation 
of two peaks is found to be less than that of a single ridge. The LED principle, however, continues to be useful 



0 0 0 0 

la: Two Peaks: TV = 4 + 2x/3 (Li), 6 (L 2 ), or 2 + 2y/3 (Loo). 



0 0 0 0 

lb: One Ridge: TV = 6 + y/Z (Li), 7 (L 2 ), or 5+3\/3 (L TO ). 


Figure 1: Breakdown of TVD: The One Ridge Case is Less Oscillatory than the Two Peaks Case. 


for multi-dimensional problems on both structured and unstructured meshes. Positivity conditions of the type 
expressed in equations (2) and (3) lead to diagonally dominant schemes, and are the key to the elimination of 
improper oscillations. The positivity conditions may be realized by the introduction of diffusive terms or by 
the use of upwind biasing in the discrete scheme. Unfortunately, they may also lead to severe restrictions on 
accuracy unless the coefficients have a complex nonlinear dependence on the solution. 

Following the pioneering work of Godunov [10], a variety of dissipative and upwind schemes designed to 
have good shock capturing properties have been developed during the past two decades [38, 6, 25, 26, 34, 31, 
11, 30, 40, 2, 15, 46, 13, 45, 5]. If the one-dimensional scalar conservation law 


dv d x _ 
dt + = 0 


( 4 ) 
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is represented by a three point scheme 


dvj . t x , 

= c i+i ^+> " v >) + c j-* to- 1 ” «>) > 

the scheme is LED if 

< + i > 0 , C ~ 1 > 0 . ( 5 ) 

J-r 2 J 2 

A conservative semidiscrete approximation to the one-dimensional conservation law can be derived by subdi- 
viding the line into cells. Then the evolution of the value vj in the j th cell is given by 


dvi 

Ax— + h- , i — h-i — 0, 
dt 3 + 2 3 2 


( 6 ) 


where /i J+ i is an estimate of the flux between cells j and j + 1 . The simplest estimate is the arithmetic average 
(fj+ 1 + fj)/ 2, but this leads to a scheme that does not satisfy the positivity conditions. To correct this, one 
may add a dissipative term and set 

h j + i = ^ (fj+i + /,) - a j + i (v j+ i -Vj). 

In order to estimate the required value of the coefficient Q^ + i, let a J+ i he a numerical estimate of the wave 
speed ^ 


du ’ 
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if Vj , -£ v , 

Vj+i-Vj J + l r 2 

if Vj +i = vj 


*1 

dv 
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and similarly 


Then 
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where 


Thus the LED condition (5) is satisfied if 
If one takes 

the diffusive flux becomes 


Av j+ i = Vj+i - Vj. 

1 | 

“i+i - 2 | a i+i ' 

1 | | 

“i+i = jP+il’ 


= I 


a > + i A «i + * 


and one obtains the first order upwind scheme 


h l= J/i ifa ;+* >0 

; + 2 1 fj + 1 if a j+i < 0 


( 7 ) 
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This is the least diffusive first order scheme which satisfies the LED condition. In this sense upwinding is a 
natural approach to the construction of non-oscillatory schemes. 

Another important requirement of discrete schemes is that they should exclude nonphysical solutions which 
do not satisfy appropriate entropy conditions [24]. These correspond to the convergence of characteristics 
towards admissible discontinuities. This places more stringent bounds on the minimum level of numerical 
viscosity [28, 41, 29, 32]. In the case that the numerical flux function is strictly convex, Aiso has recently 
proved [1] that it is sufficient that 


+ a j + i , esign(v J+ i - i>,)j 


for € > 0. Thus the numerical viscosity should be rounded out and not allowed to reach zero at a point where 
the wave speed a(u) = approaches zero. This justifies, for example, Harten’s entropy fix [11]. 

Higher order schemes can be constructed by introducing higher order diffusive terms. Unfortunately these 
have larger stencils and coefficients of varying sign which are not compatible with the conditions (2) for a 
LED scheme, and it is known that schemes which satisfy these conditions are at best first order accurate in 
the neighborhood of an extremum. It proves useful in the following development to introduce the concept of 
essentially local extremum diminishing (ELED) schemes. These are defined to be schemes which satisfy the 
condition that in the limit as the mesh width Ax ► 0, local maxima are non-increasing, and local minima are 
non-decreasing. 


2.2 High resolution switched schemes: Jameson-Schmidt-Turkel (JST) scheme 

Higher order non-oscillatory schemes can be derived by introducing anti-diffusive terms in a controlled manner. 
An early attempt to produce a high resolution scheme by this approach is the Jameson-Schmidt-Turkel (JST) 
scheme [22]. Suppose that anti-diffusive terms are introduced by subtracting neighboring differences to produce 
a third order diffusive flux 

d i+i = a >+i \ ( A ";+f + Av ;-*) } > ( g ) 

which is an approximation to The positivity condition (2) is violated by this scheme. It proves 

that it generates substantial oscillations in the vicinity of shock waves, which can be eliminated by switching 
locally to the first order scheme. The JST scheme therefore introduces blended diffusion of the form 


d i+i 


= +f + ) 4 A »>-4 


"4+i i Av i+i ~ 2Av iH + Av ;-i) ’ 


( 9 ) 


The idea is to use variable coefficients and which produce a low level of diffusion in regions where 

3 T 3 J t 2 

the solution is smooth, but prevent oscillations near discontinuities. If is constructed so that it is of order 

Ax 2 where the solution is smooth, while is of order unity, both terms in will be of order Ax 3 . 

The JST scheme has proved very effective in practice in numerous calculations of complex steady flows, 
and conditions under which it could be a total variation diminishing (TVD) scheme have been examined by 


Swanson and Turkel [39]. An alternative statement of sufficient conditions on the coefficients + , and 


for the JST scheme to be LED is as follows: 


e (4) 


Theorem 1 (Positivity of the JST scheme) 

Suppose that whenever either Vj+\ or vj ts an extremum the coefficients of the JST scheme satisfy 





= 0 


(10) 
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Then the JST scheme is local extremum diminishing (LED). 

Proof: We need only consider the rate of change of v at extremal points . Suppose that Vj is an extremum. 
Then 


<<V'5V° 


and the semi-discrete scheme (6) reduces to 

Av: + i - (V 2 \ -f Ai^_j 

di + i 2 J + V J + 2 \J~2 2 J 2 ) J 2 

and each coefficient has the required sign. □ 

In order to construct and with the desired properties define 
J 2 3 2 

R(u,v) = | 


U — V 

R+H 


if u ^ 0 or v ^ 0 
0 \i u = v = 0 


(ii) 


where q is a positive integer. Then R{u y v) — 1 if u and t* have opposite signs. Otherwise R(u , v) < 1. Now set 

Qj = R(Av j + i , Ai^_i ), Q j + i = max(Qj,Q,>i) 

and 


<s t =•,+»<&♦*. <r4=4, +1 (i-o, + i>' 

a j+4 1 


,( 4 ) 


where 


(12) 


a i+i * \ 


, /? J + iis proportional to | 

At an extremum Qj — 1, since then Ai’ J+ i and A^_i have opposite signs. Elsewhere Qj < 1 and is of order 

22 . ( 2 ) 

Ax if the solution is smooth. Thus the conditions (10) for a LED scheme are satisfied, and if q > 2, e l', is 

of order Ax 2 in smooth regions. 

An alternative formulation may be derived by noting that if Avj + | and Af^.^ are of opposite sign then 
either Vj or t>j+i is an extremum. If, on the other hand, they have the same sign, then either there is no local 
extremum, or there is a local oscillation with a maximum at either Vj or t>j+i, and a minimum at the other 
point. Suppose that e^ A is required to be zero only when At; J + | and At^-_i have opposite signs. Then if 

e^ 4 \ 0, Av- , a = <j> Av,* i where <f> > 0. Similarly if e{ 4 \ i 0, At): 3 = if>Av 4 f . i where V? > 0. Thus the 

semi-discrete scheme (6) reduces to 

Az ^ = (f+i - - (<!-* + 

+% {2A^ ; + i - (1 + *)Ar,._*} - {2A« j _ i - (1 + i>)Av j + i }. 

Since <j> > 0 and \l> > 0 it follows that the scheme is LED if 


+f - 

at every mesh point j. This is satisfied by setting 


,<«> 


A*) 


.<*) _ 
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■j+i 

= R( Av j+ 

3 , Av,! ) 

2 J 2 ' 


; ( 4 ) 

j + i 

~ 2 a -' + i( 1 

(13) 

’}+ i 

1 

> - 
- 2 
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The switches defined by the formula (11) have the disadvantage that they are active at smooth extrema. In 
order to prevent this, redefine /2(ti,i>) as 


R(u, v) = 


U — V 


max(|u| + \v \ , cAx r ) 


(14) 


where e > 0 and r is a positive power. This reduces to the previous definition if |u| + |i>| > eAx r . Now in any 
region where the solution is smooth At^+i — At or Ar^ + | — At^_i are of order Aar 2 . In fact if there is a 
smooth extremum in the neighborhood of Vj or Vj+i, a Taylor series expansion indicates that A + At) J + i 
and Avj_± are each individually of order Aar 2 , since ^ = 0 at the extremum. Choose r = Then R is of 
order Ai^. It follows that if q > 2, is of order Ax, and dj + ^ is of order Ax 2 . 

Suppose, moreover, that Vj is a maximum and that A^ + i has the opposite sign to while AVj_± 

has the opposite sign to Avj+±. Then in the scheme (13) either Qj+i = 1, or jAt^ + 3 j < fAx r and 
eAx r . Similarly either Q,_i = 1, or 

J 2 


Al):_l 
3 2 


< 


Av :+i 


< e Ai r and 


At), i 

J 2 


< cAx r . Now 


Ax 


d y j __ 

~dt ~ 


{ a j+* + ~ Q-0 ~ \ a i+i } Av i+i 

{“i-i + \ a i+^ 1 ~ Qi+O + ^ a i-i } Av i~i 


~ o “i+P 1 


Thus in any case the coefficient of At^ + i is non-negative and the coefficient of Avj_i is non-positive, while if 
Qj+± < |a^ + J < cAx r and if Qj < 1, Ia^.s 


Ai);_i > 0. It follows that 

3 i 


< cAx r . Also since v 3 is a maximum, At^+i < 0 and 


. dvj 1 . . . r 

Ax nt - i( a i+± + a :-O eAx ■ 


dt - 2' 

has the same sign as Avj_^ % then it produces a negative contribution to . In any case, therefore, 

if vj is a maximum < B , and similarly if vj is a minimum ^ > — J5, where B — ► 0 as Ax as long as 
r > 1. Thus the scheme (13) is essentially local extremum diminishing (ELED). A similar argument shows 
that the scheme (11) is also ELED provided that 


«5 a +i+ 2< i 4 + **| o i+i> 


which is the case if + ^ > \&j+^- 


2.3 Symmetric limited positive (SLIP) scheme 

An alternative route to high resolution without oscillation is to introduce flux limiters to guarantee the sat- 
isfaction of the positivity condition (2). The use of limiters dates back to the work of Boris and Book [6]. 
A particularly simple way to introduce limiters, proposed by the author in 1984 [15], is to use flux limited 
dissipation. In this scheme the third order diffusion defined by equation (8) is modified by the insertion of 
limiters which produce an equivalent three point scheme with positive coefficients. The original scheme [15] 
can be improved in the following manner so that less restrictive flux limiters are required. Let L{u,v) be a 
limited average of ti and v with the following properties: 


PI. L(u, i/) = L( v, u) 

P2. L(au, av) = oL(u, t>) 
P3. L(u, u) — u 
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P4. L(u y v) = 0 if u and v have opposite signs: otherwise L(u,v) has the same sign as u and v . 

Properties (P1-P3) are natural properties of an average. Property (P4) is needed for the construction of a 
LED or TVD scheme. 

It is convenient to introduce the notation 


(f)(r) = L(l,r) = L(r, 1). 

where according to (P4) > 0. It follows from (P2) on setting a = ^ or ^ that 

*<*.•> =*0 «=*(;)«■ 

Also it follows on setting v = 1 and u = r that 

<j>(r) = r<j> . 

Thus, if there exists r < 0 for which <f>(r) > 0, then <f> (£) < 0. The only way to ensure that <^(r) > 0 is to 
require <j>(r) = 0 for all r < 0, corresponding to property (P4). 

Now one defines the diffusive flux for a scalar conservation law as 



d j+± = a f+i { At 5+* 


(15) 

Also define 

As. ,3 

_+ _ 3 + 3 

1 

1 

1 

> 

1 

K)lW 



r~ — — r 

J 3 

Then, the scalar scheme (6) reduces to 


1 A 1 A 

= ~2 a iH Av j + i ~ 2 a J-i Av J-i 
+« ; + i (a^ ; + i - <i>(r+)Av j _ i ) 

~ a j-i ( Av s-i 



a i+i ~ \ a i+i + )} 

Q i-^ + + «i+^( r+ )} Av i-i 


(16) 


Thus the scheme satisfies the LED condition if a j+ i > ~ |a ;+ i| for al1 J. and <M r ) > 0. which is assured 
by property (P4) on L. At the same time it follows from property (P3) that the first order diffusive flux is 
canceled when Av is smoothly varying and of constant sign. Schemes constructed by this formulation will be 
referred to as symmetric limited positive (SLIP) schemes. A variation is to include the coefficient a ; + i in the 
limited average by setting 


d j+$ ~ ( 17 ) 

-L (a j+ iAv j+ z,a j _iAv j _^J . 

It is easily verified that the argument remains valid. These results may be summarized as 
Theorem 2 (Positivity of the SLIP scheme) 

Suppose that the discrete conservation law (6) contains a limited diffusive flux as defined by equations (15) or 
(17). Then the positivity condition (7), together with the properties (P1-P4) for limited averages , are sufficient 
to ensure satisfaction of the LED principle that a local maximum cannot increase and a local minimum cannot 
decrease. □ 
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The construction benefits from the fact that the terms involving <f>(r~) and <£(r+) reinforce the positivity 
of the coefficients whenever <f> is positive. Thus the only major restriction on L(u,v) is that it must be zero 
when it and v have opposite signs, or that <t>{r) — 0 when r < 0. If Avj+i and Ai^_i have opposite signs then 
there is an extremum at either j or j -f 1. In the case of an odd-even mode, however, they have the same sign, 
which is opposite to that of At; J + ^, so that they reinforce the damping in the same way that a simple central 
fourth difference formula would. At the crest of a shock, if the upstream flow is constant then A Vj_i = 0, 
and thus At> ; - + | is prevented from canceling any part of At>j+^ because it is limited by Avj_± . 

A variety of limiters may be defined which meet the requirements of properties (P1-P4). Define 

S(u, v ) = \ { si « n ( u ) + s '8 n ( v )} 


so that 


S{u,v) 


1 if u > 0 and v > 0 
< 0 if u and v have opposite sign 

— 1 if it < 0 and v < 0 


Three limiters which are appropriate are the following well-known schemes: 


1. Minmod: 

L(u, v) = 5(u, v ) min(|u|, \v\) 


2. Van Leer: 


L(u, d) = S(u, v) 


2MM 
M + M 


3. Superbee: 


L(u , u) = S(u , v) max {min (2|u|, \v \) , min (|u|, 2|v|)} 


Another formulation is simply to limit the arithmetic mean by some multiple of the smaller of |u| and 


4. a- mean: 


L(u, u) = S(u, v) min 


| u + ^1 
2 



With the present construction the first three of these limiters are unnecessarily stringent. Superbee, for 
example, could be relaxed to 
o-bee: 

L(u y v) = 5(u, v) max {min (a|u|, |v|) , min (|u|, Qf|v|)} 

which reduces to minmod when a = 1, and is less stringent then Superbee when a > 2. Superbee differs 
from the other limiters in that it introduces a larger amount of antidiffusion than that needed to cancel the 
diffusion, + f), when | < J < 3. The resulting negative diffusion tends to produce artificial compression 
of discontinuities. 

In order to produce a family of limiters which contains the first two limiters it is convenient to set 


L(u,v) 


1 

2 


D(u, a)(u + v) 


where D(u, t>) is a factor which should deflate the arithmetic average, and become zero if u and v have opposite 
signs. Take 


D(u } u) = 1 — R(u, v ) — 1 — 


u — v 


19 


M + M 


(18) 


where is the same function that was introduced in the JST scheme, and q is a positive integer. Then 

D(u y v) = 0 if u and v have opposite signs. Also if q = 1, L(u,v) reduces to minmod, while if q = 2, L(u,v) 
is equivalent to Van Leer’s limiter. By increasing q one can generate a sequence of limited averages which 
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approach a limit defined by the arithmetic mean truncated to zero when u and v have opposite signs. When 
the ratio r = ^ is extreme the limiter approaches the asymptotic value 

<j>(r) = L(1 , r) — ► q as r — ► oo. 

When the terms are regrouped it is apparent, moreover, that the JST scheme with the second definition of the 
switch, equation (13), and the SLIP scheme with this limiter are identical. This unifies the two formulations. 

As in the case of the JST scheme, the SLIP scheme can be relaxed to give an essentially local extremum 
diminishing (ELED) scheme which is second order accurate at smooth extrema by introducing a limited average 
with a threshold through the definition 


v ' I max(|u| 4- |v| , e Ax r ) | 

where r = q > 2. The effect of this “soft limiter” is not only to improve the accuracy: the introduction of a 
threshold below which extrema of small amplitude are accepted also usually results in a faster rate of conver- 
gence to a steady state, and decreases the likelyhood of limit cycles in which the limiter interacts unfavorably 
with the corrections produced by the updating scheme. In a scheme recently proposed by Venkatakrishnan a 
threshold is introduced precisely for this purpose [43]. 

2.4 SLIP schemes on multi-dimensional unstructured meshes 


2 



Figure 2: Cell Surrounding Vertex o. 


Consider the discretization of the scalar conservation law (1) by a scheme in which v is represented at the 
vertices of a triangular mesh, as sketched in Figure 2. In a finite volume approximation (1) is written in 
integral form as 

J vdS + j> (f(v) dy - g(v) dx ) = 0 , 

and this is approximated by trapezoidal integration around a polygon consisting of the triangles with a common 
vertex, o, say. 

Thus (1) is discretized as 

s^- + \ 51 {(/* + A-0 (j tk - yk-\) - ( 9k + 9 k-i)(x k - = o 

k 

where /* = /( v k ), g k = g{v i fc ), 5 is the area of the polygon, and k ranges over its vertices. This may be 
rearranged as 

5 ^ + E(A Ay* - = 0 
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where 


Ax* = - (xjt+i - x*_i) , Ay* = - (y*+i - y*-i) . 

Following, for example, References [16] and [21], this may now be reduced to a sum of differences over the 
edges ko by noting that £2* Ax* = Ay* = 0. Consequently f„ and g„ may be added to give 

S ^t + _ U)^yk - (gk - y 0 ) Ax*} = 0. (19) 

Define the coefficient ajt 0 as 


0>ko ~ 


< 



V* ^ Vo 


V* = v 0 


and 


Av ko = v k - tv 

Then equation (19) reduces to 

s ~^f + XI ak ° Avk ° = 

k 

To produce a scheme satisfying the sign condition (2), add a dissipative term on the right hand side of the 
form 

]Pa* 0 Au* o , (20) 

* 


where the coefficients a ko satisfy the condition 


&k o ^ |^fco | • 


(21) 


These simple schemes are far too dissipative. Antidiffusive terms may be added without violating the 
positivity condition (2) by the following generalization of the one-dimensional scheme. Considering again the 
scalar case, let l ko be the vector connecting the edge ko and define the neighboring differences 

A + t >* 0 = Iko-V+v, A'Vko = Iko-V-v, 

where are the gradients of v evaluated in the triangles out of which and into which l ko points, as sketched 
in Figure 3. Arminjon and Dervieux have used a similar definition [3]. 



Figure 3: Edge ko and Adjacent Triangles. 


It may now be verified that 


A+V ko = e pJt (V p - V k ) + t qk (V q - V k ) 
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and 


A V}e 0 — i or (l> 0 V r ) -f- t 0 s (^o ^s) i 

where the coefficients c p *, e q k> <-or and c os are all non-negative. Now define the diffusive term for the edge ko 
as 

dfco — &ko { AtJk 0 E (A"^t>k0, A Vko) } j (^) 

where L(u,v) is a limited average with the properties (P1-P4) that were defined in Section 2.3. In considering 
the sum of the terms at the vertex o write 

L(A + v ko , = <j>(r£ 0 ) A~v ko , 


where 

. A^v^q 
r *°~ 

Then, since the coefficients c or and e os are non-negative, and 4>(r% 0 ) is non-negative, the limited antidiffusive 
term in (22) produces a contribution from every edge which reinforces the positivity condition (2). Similarly, 
in considering the sum of the terms at k one writes 

L(A + v k0 , A~v ko ) = if>(r~ o )A + v ko , 


where 


r 


ko 


A Vko 

A+t>jbo 7 


and again the discrete equation receives a contribution with the right sign. One may therefore deduce the 
following result: 


Theorem 3 (Positivity Theorem for Unstructured Meshes) 

Suppose that the discrete conservation law (19) is augmented by flux limited dissipation following equations 
(20) and (22). Then the positivity condition (21), together with the properties (P1-P4) for limited averages , 
are sufficient to ensure the LED property at every interior mesh point. □ 

Note also that if this construction is applied to any linear function v then 


Av ko - A+Vko = A Vkoy 


with the consequence that the contribution of the diffusive terms is exactly zero. In the case of a smoothly 
varying function v y suppose that l ko -Vv ^ 0 and the limiter is smooth in the neighborhood of rf 0 = 1. Then 
substitution of a Taylor series expansion indicates that the magnitude of the diffusive flux will be of second 
order. At an extremum the antidiffusive term is cut off and the diffusive flux is of first order. 


2.5 Upstream limited positive (USLIP) schemes 

By adding the anti-diffusive correction purely from the upstream side one may derive a family of upstream 
limited positive (USLIP) schemes. Corresponding to the original SLIP scheme defined by equation (15), a 
USLIP scheme is obtained by setting 


if Oj + i > 0, or 




d i+± 


= «;+* 

= “i+i { A *;+± 


— L ( A r J+ ^ , Avj_fj } 

- L ( Av i+^ Av i+i)} 
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ifa J + i <0. If « j+ i = 2 


one recovers a standard high resolution upwind scheme in semi-discrete form. 


Consider the case that a. , i >0 and a,_i >0. If one sets 

3 -r-j 3 2 


r+ _ Av i+i r -_ At, i-l 

r — — , r = 


A v** 


Av ,-i ' 


the scheme reduces to 

Ax^ = {«i(r + )a J+ i + (2 

To assure the correct sign to satisfy the LED criterion the flux limiter must now satisfy the additional 
constraint that <j>(r) < 2. 

The USLIP construction can also be implemented on an unstructured mesh by taking 


dko = |Oiko| {At?*,, - L (Aut 0 , A Vjto)} 


if ajfc 0 > 0 and 

dko = |a*o| {At>* 0 - L (Avk„,A + Vk 0 )} 

if dko < 0- Let and denote sums over the edges meeting at the vertex o for which dko > 0 and 

A + v ko 


dko < 0. Define 


+ _ Avko 

Tko ~ A -v ko ' Tko ~ 


A Vko 


Then 


,dv 0 
' dt 


= - a ko<t> (rto) a Vko 

- °ko ( 2 ~d> (r-ko)) Avfco 


and substituting the formula for A Vk 0 the coefficient of every difference A t>* 0 is found to be nonnegative, 
with the consequence that the scheme is LED. 


2,6 General construction of higher order SLIP schemes 

Schemes of any desired order of accuracy in regions where the solution does not contain extrema can be 
constructed by the following general procedure. Suppose that the scalar conservation law (1) is approximated 
in semi-discrete form by the low and high order schemes 


and 






(23) 


(24) 


where the low order scheme has positive coefficients and is local extremum diminishing (LED). Define an 
anti-diffusive flux as 

A j + ±=fH 1 + l-h, +i ' 

and in order to define a limited corrective flux fc j+ i let + \ be a bound determined from the local slopes as 


J i + { 


minmod (Av J+ | , Av j+ x, Av 7 _i) | , 


where minmod (u, t>, u;) = 0 if u, v , and w do not have the same sign, and otherwise 

minmod (u, v,w) = S min (|ti|, |v|, \w \) , 
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where S is the sign of u, v , and w. Set 

fo i+ 1 = s '8 n (^+i) min (|^j+i|.^+i B i+i) - 
where , i > 0. The SLIP scheme is now defined as 

J ' 2 

dv ■ 

Ax — + hi i i — h-> = 0, 

dt 1 + i 3 * 

where 


(25) 


(26) 


(27) 


h i+\ = h,^ +fc, +i - 

Thus, it reduces to the high order scheme when the limiters are not active. It is important that the limiter 
depends only on the magnitude of the local slopes, and not their sign, so that the correction can have either 


sign. 


In order to prove that the general SLIP scheme is LED, note that the low order LED scheme can be written 


Ax 


= ct i At),- . 1 — C l Al): 1 , 
dt J+i J + 2 J-i 3 2 


where ct . > 0 and c. ± >0. Now, 

J + 2 J 2 

min(|A_, + i 0 <-y J+ x < 0 }+ i- 


Also, 


where 


and 


*j = ^+i A ^+i =^+iA v i-i> 


*,+4 =sign(At> i+ i) 


o < <AJ + a < 1, o < <t> j+i < 1, 0 < <t>- + ± < 1 

since B , • , ± = 0 if Ad, , 2 , Ad, , 2 , and Ad- 1 do not all have the same sign. Define 

JT 2 J T 2 JT 2 J 2 

^,+4 = si 6 n ( A i+i)- 

Then, since + = Bj_i = 0 unless = s ; - _ ^ , 


~fc 3+i = 


1 

+ 2 

1 

~2 

1 

+ 2 
_1 

“2 




S j-| S V* 




s >+i 5 *, + *| 

so that the correction reinforces the positivity of ct , , and c7 This provides the proof of 

j + 2 J 2 


Theorem 4 (Positivity of the General SLIP Scheme) 

The semi- discrete scheme defined by equations (25-27) is LED if the low order scheme (23) is LED . □ 

The key idea in this proof is that the correction Aj+i may be associated with either Ad^_i or A ^ + i depending 
on whether it has the same or the opposite sign as Ad^.i and Ad j + ^. 

The idea of blending high and low order schemes to produce a limited anti-diffusive correction is similar to 
that used in Zalesak’s generalization of flux corrected transport (FCT) [47]. With FCT the anti-diffusion is 
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introduced in a separate corrector stage, whereas in the present scheme it is integrated in the construction of 
the numerical flux. This brings it within the framework of a general theory of LED schemes, and facilitates 
its extension to treat systems of equations by the introduction of flux splitting procedures. 

To prevent the scheme collapsing to the low order scheme at smooth extrema, modify the definition of 

j-r 2 

to 

Bj + i = max j|minmod + i , Ai^ + ^ , Ai^_±^ j ,eAx r | . 

Now either 

s j+i B j + \ = *i+i Av i+l = 4> j+ iAv j+i = 


or 






< P j+ $t Ax T 

Suppose that the low order scheme is essentially local extremum diminishing (ELED). Then if Vj is a maximum 


fh-i < P > 1 

It follows that the corrected scheme satisfies 

^ < KlAx”-' + (0 j+i + 0 j _±)tAx r ~ l 
Taking r > p, there is a constant K such that: 


< KAx p 
dt 


Similarly if vj is a minimum 


dvj 

~dt 


> -KAx p 


Therefore the corrected scheme is ELED. Also if the low order scheme is of order L, then in a region where 
the flow is smooth 


A i+i 


0(Ax l ) 


since the leading error term is cancelled by the high order scheme. If r > 1, and L > r, then when the mesh is 
sufficiently fine + ^ will not be limited by B near a smooth extrema, so the accuracy of the high order 
scheme will be recovered. 

As an example of the general SLIP construction suppose that the numerical flux has the form 


~ 2 ^ +1 + 

where for the low order scheme 


d :+i 


~ a i+i Av i+* 


a iH Z o 




and for the high order scheme 


d i+i = ( At; >+* ~ ^ Av i+i ~ \ Av i-l) ■ 

These are just the diffusive fluxes which are used in the switched JST scheme described in Section 2.2. The 
anti-diffusive flux in the SLIP scheme is now 



In this case the bound B j+ ^ need only depend on the smaller of | Av ;+ j | and |At^_i | , provided that A^ + i 
and Av ; _^ have the same sign, leading to the first SLIP scheme with a-mean as the limiter. Here the SLIP 
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construction provides an alternative switching procedure to the sensor in the JST scheme, such that the LED 
property is enforced. 

To construct a sequence of successively higher order SLIP schemes one may start by constructing a second 
order scheme SLIP 2 , say by taking 


f (1) 

JL 




fO) 

JH 


+ i 


= ,(/j+i +/j)-a J+ iAt> j+i 


2 (/i+i + fj) 


* 15 . 


( Av j+f _2AtJ ;+i + 


(i) _ 


A 1 ) A 1 ) 


- * 0 ) 


- 

— J L 


+ i 




Then one may repeat the procedure, taking 



f< 2 > 

kl >+i 





f (2) 

2 (/j+1 A- /;) “ 




1 , 

6“ J+ i ' 


~ 4Av ,+ i 

H-6Av J + i -4Av } _i + At^aJ 

where 









A /; + i = 

A U 

+.-A fj, 

and 







A w _ 
A *+i ~ 

A 2) 

1 

- f (2) 




= 

A2) 

+ sign ^.4 

(2) ^ 



The resulting scheme, which may be conveniently labelled SLIP4, is fourth order accurate when the limiters are 
inactive. The procedure may then be iterated. The correction is of order Ar 2 , and subsequent 

corrections are of correspondingly higher order. Thus they are progressively less likely to be limited by the 
bound R , 1 . 

J-r 2 


2.7 General SLIP scheme on unstructured meshes 

The general SLIP construction may also be implemented on unstructured meshes. With the notation of 
Figure 4, let fi ko and fn ko be low and high order fluxes along the edge ko. Define the anti-diffusive flux along 
this edge as 

A ko = fH ko - h ko ( 28 ) 

and the limited corrective flux as 

f Cko = sign(^jbo) min(|>Uo| , PkoBko ) , (29) 

where f3 ko > 0 and Bko is a bound determined by the local slopes. In order to define Bko let / and n be any 
vertices neighboring o and k such that 

sign(Ai> n fc) = sign(Ai; / bo), sign(Av 0 j) = sign(Aujt 0 ). 

If there is no such vertex n then A: is a local extremum , and if there is no such vertex / then o is a local 
extremum. In either case set Bk 0 = 0. Otherwise set 

B ko = min(|Av„jk| , \Av ko \ , |Au 0 /|) • (30) 
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I 


Figure 4: Edge ko and Adjacent Edges. 


The flux along the edge ko for the SLIP scheme is now defined as 

fko = h ko + /c fco * 

It may be verified that the scheme can be expressed in terms of differences between the vertex o and its 
nearest neighbors with non-negative coefficients by adapting the one-dimensional derivation of the last section 
in the same way that the one-dimensional derivation of Section 2.3 was adapted to the unstructured mesh in 
Section 2.4. This result may be stated as 


Theorem 5 (Positivity of the General SLIP Scheme on Unstructured Meshes) 

If the discrete conservation law (19) is augmented by the diffusive flux fi ko and fc ko defined by equations (28- 
SO), then the scheme is LED at every interior point. □ 

The construction requires the identification of any three edges lo } ok and kn along which the solution is 
monotonically increasing or decreasing. If Avk 0 > 0 one could search for vertices / and n which maximize 
At> n * and A v c i, but since the number of edges meeting at a given vertex can be very large, this procedure 
could be expensive, and one might prefer to apply the test to the edges nk and ol most nearly aligned with 
the edge Jbo. 

2.8 Fully discrete LED schemes 

When a discrete time stepping scheme is introduced to produce a fully discrete scheme, let a superscript n 
denote the time level, and suppose that 

k 

A Taylor series expansion now shows that if the discrete scheme corresponds to a differential equation with no 
source term, then 

E 6 ** = 1. 

k 

Also 

b n+1 | ^ Ei 6 **i m * ax i v *l 

k 

Therefore the solution has a nonincreasing norm if 

Em<i- 

k 
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Upper Threshold for j+1 



Figure 5: Thresholds for a Fully Discrete LED Scheme. 


These two conditions can be satisfied only if b jk > 0 for all j, k. The LED principle may now be expressed 
by requiring that < t> ; n if vj is a maximum and u" +1 > vj if vj is a minimum, while if v? is not an 

extremum, must be within thresholds defined by the maximum and minimum values of v k at the nearest 
neighboring points. This will be the case if both bjk > 0 and bjk — 0 if j and k are not neighbors. The 
following result is immediate. 


Theorem 6 (Positivity of Fully Discrete Schemes) 

Suppose that the semi-discrete scheme 

is LED. Then a time step At > 0 can be found such that the corresponding forward Euler scheme 

V ? +1 = + At *(«*-«") 

is LED. 

Proof: The coefficients of the discrete scheme are 

bj k — At dj k , k ^ j 
bjj = 1 — At ^ ajk 

The off diagonal coefficients bjk inherit the positivity of the coefficients of the original scheme, while bjj >0 if 


□ 

Given a forward Euler scheme that satisfies the positivity conditions, Shu has devised a procedure for con- 
structing higher order multi-stage time stepping schemes which preserve these conditions under an appropriate 
restriction of the time step [35]. Second and third order schemes can be constructed without any modification 
of the space discretization. 
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In the case of the one dimensional conservation law a forward Euler scheme has the form 


„n+ 1 


At 




and one must consider the interaction of the flux /i J + i with the solution at the points j + 1 and j. Both the 
SLIP and USLIP constructions remain valid under a constraint on the time step which depends on the choice 
of the limiter. Suppose that the flux limiting function <j>(r) is bounded, <j>(r) < </> max . Then the limit on At 
becomes smaller as <^ ma x is increased. 

In the case of the SLIP scheme, for example, equation (16) is replaced by 


•r = »” 


+ 


Af 1 

r 1 , a 

| 

Ax 1 

l a j+i~ 2 a i+i +a i-i^ r 

f K+i - «") 

At J 

r 1 , .1 

1 

Ax 1 

( 0 J-i + 2 a i-i + a f+i^( r )j 

f K - v j- 1 ) 


The coefficients of i >” +1 and i?"_j are non-negative if <*j+i > \ | a j + £ • If one takes the minimum level of 
diffusion a, . i = 77 

j » 2 i 


j+\ 


the coefficient of v” is 


L _ , A t 
bjj - 1 2Ax 


{| a i + i I ~ a j +i + ^( r ) | a >-i | + | a i-§ | + a i-i + | a j+i } 


Consider the case when both a ; + x and are positive. Take a = max(a J + x, a^i). Then bjj > 0 if 

or in the case of the limiter (18), a 


3 Systems of conservation laws 

3.1 Flux splitting 

Steger and Warming [38] first showed how to generalize the concept of upwinding to the system of conservation 
laws 

dw 0 

_ + _ /H = 0 (31) 

by the concept of flux splitting. Suppose that the flux is split as / = / + +/~ where and have positive 
and negative eigenvalues. Then the first order upwind scheme is produced by taking the numerical flux to be 

^i+i = ft + fj+ 1 - 

This can be expressed in viscosity form as 


l j+\ 


- +2 (/i+i + ft) 2 (ft+i - ft) 

+ o (fj+ 1 + fj) + 9 (//+ 1 “ /f) 


where the diffusive flux is 


1 


d j+ i = -A(f+-f-) jH . (32) 

Roe derived the alternative formulation of flux difference splitting [34] by distributing the corrections due to 
the flux difference in each interval upwind and downwind to obtain 

+ (/>+» - fiT + (fi - fj- 0 + = 0 , 

where now the flux difference /j+i — fj is split. The corresponding diffusive flux is 


^J = iK +j -A/- + ,) 
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Following Roe’s derivation, let Aj + i be a mean value Jacobian matrix exactly satisfying the condition 


/j+i — ^ + -1(^ + 1 — w j)- 

Then a splitting according to characteristic fields is obtained by decomposing Aj + ± as 


Aj H =ThT~\ 


where the columns of T are the eigenvectors of and A is a diagonal matrix of the eigenvalues. Then 

A =TA ± r 1 At» i+l . 

Now the corresponding diffusive flux is 


n A j+1 ( W j + 1 - W })> 


where 


A i+h 


T|A| T~ l 


and | A | is the diagonal matrix containing the absolute values of the eigenvalues. 
Simple stable schemes can be produced by the splitting 


(/;+ 1 — />) = 2^+ 1 ~ 1 -Wj), 

which satisfies the positivity condition on the eigenvalues if a J+ i > ^ max |a(j4j + i) 
scalar diffusive flux 

d i+\ = a ;+i Aw > + |- 


and corresponds to the 


(33) 


Characteristic splitting has the advantage that it allows a discrete shock structure with a single interior point. 
The simple scalar diffusive flux (33) is computationally inexpensive, and combined with the high resolution 
switched scheme captures shock waves with about three interior points. 


3.2 Construction of convective upwind and split pressure (CUSP) schemes 


Discrete schemes should be designed to provide high accuracy in smooth regions in combination with oscillation- 
free shocks at the lowest possible computational cost. This in turn requires both economy in the formulation, 
and in the case of steady state calculations, a rapidly convergent iterative scheme. The convective upwind 
and split pressure (CUSP) scheme described below meets these requirements, while providing excellent shock 
resolution at high Mach numbers. When very sharp resolution of weak shocks is required, the results can be 
improved by characteristic splitting with matrix diffusion using Roe averaging. 

Consider the one-dimensional equations for gas dynamics. In this case the solution and flux vectors appearing 
in equation (31) are 


w = 


P ^ 
pu 

pE ) 


f = 


pu \ 

pu 2 + p , 

puH J 


where p is the density, u is the velocity, E is the total energy, p is the pressure, and H is the stagnation 
enthalpy. If 7 is the ratio of specific heats and c is the speed of sound 


p = (7 - 1 )p 



c 


2 


7 P 
P 
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In a steady flow H is constant. This remains true for the discrete scheme only if the diffusion is constructed 
so that it is compatible with this condition. 

The eigenvalues of the Jacobian matrix A = are u, u + c, and u — c. If u > 0 and the flow is locally 
supersonic (M = ~ > 1), all the eigenvalues are positive, and simple upwinding is thus a natural choice for 
diffusion in supersonic flow. It is convenient to consider the convective and pressure fluxes 




(o' \ 

fc = U 

pu 

= UW C , fp=\ P 1 


P" ) 

\ o ) 


separately. Upwinding of the convective flux is achieved by 


c >+4 


l i+4 


Au S +i = \M\ c i+i Auv +i , 


where M is the local Mach number attributed to the interval. Upwinding of the pressure is achieved by 


d p, + ± = sign(M) 


0 'I 


Full upwinding of both f c and f p is incompatible with stability in subsonic flow, since pressure waves with 
the speed u — c would be traveling backwards, and the discrete scheme would not have a proper zone of 
dependence. Since the eigenvalues of -|^ are u, u and yu , while those of are 0, 0 and —(7 — l)u, a split 
with 

f + = fc r = f P 

leads to a stable scheme, used by Denton [ 8 ], in which downwind differencing is used for the presure. 

This scheme does not reflect the true zone of dependence in supersonic flow. Thus one may seek a scheme 
with 



fi(M)c j+ xAw Cj+i 
MM) ( A Pj+i j , 


where f\(M) and /^(M) are blending functions with the asymptotic behavior /i(M) — » |M| and / 2 (M) — ► 
sign(M) for \M\ > 1. Also the convective diffusion should remain positive when M = 0, while the pressure 
diffusion must be antisymmetric with respect to Af . A simple choice is to take f\(M) = |Af| and / 2 (M) = 
sign(Af) for | M \ > 1, and to introduce blending polynomials in M for \M\ < 1 which merge smoothly into the 
supersonic segments. A quartic formula 


/ 1 (M) = a 0 + a 2 M 2 + a 4 M 4 , \M\ < 1 

preserves continuity of f\ and ^ at \M\ = 1 if 

_ 3 _ 1 

— 2 — <Jo 2* 

Then a c controls the diffusion at M = 0. For transonic flow calculations a good choice is a Q = while for 
very high speed flows it may be increased to A suitable blending formula for the pressure diffusion is 

/ 2 (A/) = Im(3-M 2 ), |M|<1. 

The diffusion corresponding to the convective terms is identical to the scalar diffusion of Jameson, Schmidt 
and Turkel [22], with a modification of the scaling, while the pressure term is the minimum modification 
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needed to produce perfect upwinding in the supersonic zone. The scheme retains the property of the original 
scheme that it is compatible with constant stagnation enthalpy in steady flow. If one derives the viscosity 
corresponding to the flux splitting recently proposed by Liou and Steffen [27], following equation (32), one 
finds that their scheme produces first order diffusion with a similar general form, and the present scheme may 
thus be regarded as a construction of artificial viscosity approximately equivalent to Liou-Steffen splitting. 


3.3 Multi-dimensional systems 

Schemes for structured meshes are conveniently constructed treating each mesh direction separately in a 
manner similar to the one-dimensional case. For unstructured meshes, the three-dimensional conservation law 

m + + s/ w + ai h(v) = 0 (34) 

can be treated in a manner similar to the scalar case by first expressing the convective flux balance as a sum of 

differences along edges. Consider the set of tetrahedrons containing a common edge. Then one may associate 
with that edge a vector area S which is one-third the sum of the areas of the set of faces which form one of two 
opposing umbrellas around the edge. With a notation similar to that of Figure 2 the convective flux balance 
corresponding to equation (34) at an interior mesh point may be written as 

V^ + ^( F it- F o )- s io = 0 , (35) 

k 

where the columns of F are the flux vectors /, g and h, and V is the volume of the polyhedron formed by 
the union of all the tetrahedrons with the common vertex o. Here F 0 may be added or subtracted since 

^jfcSjto = 0. Diffusion may now be added along the edges in exactly the same way as before. When the 

convective flux balance is evaluated, it is more convenient to use the sum ^2 k (F k + F 0 )'Sjt o , so that the 
convective flux along each edge needs to be calculated only once in a loop over the edges and appropriately 
accumulated at nodes k and o. 

The SLIP scheme can now be formulated with the aid of Roe’s construction [34]. Let A ko be a matrix such 
that 

A ko (^ib W 0 ) — (Ffc F 0 ) *S ko . 

Suppose that A ko is decomposed as TAT" 1 where the columns tj of T are the eigenvectors of A ko • Then 
the difference Aw = w k - w 0 is expressed as a sum ajtj of the eigenvectors, where the coefficients aj = 

(T~ l Aw). represent the characteristic variables, and the diffusive term along the edge ko is constructed as 

\Ako\ Aw = T |A| T~ l Aw. 

In order to construct a higher order scheme, an anti-diffusive flux may then be calculated by applying the 
limited averaging procedure as in equation (22) to each characteristic variable separately. 

At boundary points equations (19) or (35) need to be augmented by additional fluxes through the boundary 
edges or faces. The first order diffusive flux ot ko Av ko may be offset by subtracting an antidiffusive flux evaluated 
from the interior, taking a limited average with Av ko . 


4 Convergence acceleration for steady state calculations 

4.1 Time stepping schemes 

The discretization of the spatial derivatives reduces the partial differential equation to a semi-discrete equation 
which may be written in the form 


^ = °- 


( 36 ) 



where w is the vector of flow variables at the mesh points, and R(w) is the vector of the residuals, consisting 
of the flux balances augmented by the diffusive terms. In the case of a steady state calculation the details of 
the transient solution are immaterial, and the time stepping scheme may be designed solely to maximize the 
rate of convergence. 

If an explicit scheme is used, the permissible time step for stability may be so small that a very large number 
of time steps are needed to reach a steady state. This can be alleviated by using time steps of varying size 
in different locations, which are adjusted so that they are always close to the local stability limit. If the 
mesh interval increases with the distance from the body, the time step will also increase, producing an effect 
comparable to that of an increasing wave speed. Convergence to a steady state can be further accelerated 
by the use of a multigrid procedure of the type described below. With the aid of these measures explicit 
multistage schemes have proved extremely effective. 

If one reduces the linear model problem corresponding to (36) to an ordinary differential equation by 
substituting a Fourier mode w = e ,px ^, the resulting Fourier symbol has an imaginary part proportional to 
the wave speed, and a negative real part proportional to the diffusion. Thus the time stepping scheme should 
have a stability region which contains a substantial interval of the negative real axis, as well as an interval 
along the imaginary axis. To achieve this it pays to treat the convective and dissipative terms in a distinct 
fashion. Thus the residual is split as 


R(w) - Q{w) + D(w), 

where Q(w) is the convective part and D(w) the dissipative part. Denote the time level nAt by a superscript 
n. Then the multistage time stepping scheme is formulated as 

^(n + 1,0) _ u ,n 

u>< n+1 '*> = w n -a k At (q(*- 1 ) + D<*- 1 >) 

U, n + 1 = W ( n + 1 ' m ), 


where the superscript k denotes the k - th stage, a m = 1, and 


= Q(w n ), D^ = D{w n ) 

Q (k) = Q (u/ n+1 '*>) 

D (k] = 0 k D + (1 - p t )D (k - 1) . 

The coefficients c** are chosen to maximize the stability interval along the imaginary axis, and the coefficients 
/?* are chosen to increase the stability interval along the negative real axis. 

These schemes do not fall within the standard framework of Runge-Kutta schemes, and they have much 
larger stability regions. Two schemes which have been found to be particularly effective are tabulated below. 
The first is a four-stage scheme with two evaluations of dissipation. Its coefficients are 


<*i = 

1 

3 

0i 

C*2 = 

4 

15 
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5 
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<*4 = 

1 
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( 37 ) 
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The second is a five-stage scheme with three evaluations of dissipation. Its coefficients are 


\ 

0i = 1 

«2 = i 

02=0 

«3 = f 

0 3 = 0.56 

— KN 

II 

TT 

S 

04=0 

a 5 = 1 

05 = 0.44 


(38) 


4.2 Multigrid 

The multigrid scheme is a full approximation scheme defined as follows [14, 16]. Denote the grids by a subscript 
Jb. Start with a time step on the finest grid k = 1. Transfer the solution from a given grid to a coarser grid by 
a transfer operator Pk,k-u so that the initial state on grid k is 

(0) D 

= /Tfc.fc-lU't-l- 

Then on grid k the multistage time stepping scheme is reformulated as 

w[ q + l) = u4 0) - a n At + Gkj , 

where the residual R^ is evaluated from current and previous values as above, and the forcing function G* 
is defined as the difference between the aggregated residuals transferred from grid k — 1 and the residual 
recalculated on grid k. Thus 

G k = Q k ,k-iR(w k - 1 )-R(w[ 0) ), 

where Qk,k-\ is another transfer operator. On the first stage the forcing term G* simply replaces the coarse 
grid residual by the aggregated fine grid residuals. The accumulated correction on a coarser grid is transferred 
to the next higher grid by an interpolation operator h-\,k so that the solution on grid Ar — 1 is updated by 
the formula 

U>k - 'll = “’*-1 + Ik-1, k (w k ~ t4 0) ) • 

The whole set of grids is traversed in a W -cycle in which time steps are only performed when moving down 
the cycle. First order numerical diffusion is always used on the coarse grids, and in cases when characteristic 
splitting is used on the fine grid, simple scalar diffusion is used on the coarse grids. 

5 Numerical Results 

Extensive numerical tests have been performed with schemes based on the theory of Sections 1-4. Some 
results are presented here to illustrate the performance in practice of both the CUSP scheme and schemes 
using characteristic splitting. 

5.1 Shock tube 

In order to verify that the general higher order SLIP construction is effective in preventing oscillations in 
unsteady flow, Figure 6 presents a calculation of the Sod problem [36] for a shock tube. This calculation 
was performed with characteristic splitting, using the fourth order SLIP scheme formulated in Section 2.6. 
In the one-dimensional case single stage time stepping schemes of second or higher order similar to the Lax- 
Wendroff scheme can be derived at little cost in complexity by the successive substitution of space derivatives 
for time derivatives. Here, since the purpose was to verify the LED property of the SLIP scheme, a simple 
forward Euler time stepping scheme was used with a time step corresponding to a Courant number of 1. The 
computed results are superposed on the exact solution. The shock wave and expansion are very well resolved. 
The contact discontinuity is less sharply resolved, as is to be expected because of the absence of a natural 
compressive effect at a contact discontinuity. 
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5.2 Steady flow calculations using the CUSP scheme 


A variety of multigrid calculations using the CUSP scheme are presented in Figures 7-11, including both 
transonic and hypersonic flows. The scheme was implemented with a JST switch which was formulated with 
pressure gradients. Using subscripts i f j to label the mesh cells, let j be the diffusive flux calculated by 
the CUSP scheme from the states and Wij . Then the final diffusion is 


where 


II 

■*-5 

-K'* 

+ 

c (2) 
'*+ 2’- 

( 4 ) 

“ € , + l j( e *‘+|j ~ 2e * + i,j + 1 

e (2) 


aR i+$,j’ e !+ij = max (°- 1 - 


= 

max( , Ri+ij , ftij , Rj—ij) 

Rij 


RIApi+ij.Api.ij) 

A P«'+i,i 

- 

Pi+ij ~ Pij 


and R(u, v) is the function defined by equation (14) with q = 3, r = |. In hypersonic flow a = 1, corresponding 
to pure upwinding, but in transonic flow the shock resolution is improved by taking a = .625. In both cases 
3 is taken as 2. 

With this construction the role of the high order diffusion is to provide global damping of oscillatory modes 
which would otherwise inhibit convergence to a steady state, while the role of the first order diffusion is to 
control oscillations near discontinuities. Numerical experiments with multigrid acceleration confirm that the 
rate of convergence to a steady state is essentially the same when the first order diffusion is eliminated, but 
large pre- and post-shock oscillations appear in the solution. On the other hand the multigrid scheme will not 
converge if the global diffusion is eliminated. 

Figures 7-9 show transonic solutions for three different airfoils, calculated on 160x32 meshes with O- 
topology, and each of which is essentially converged in 12 multigrid cycles. The five stage time stepping 
scheme (38) was used, and the work in each cycle is about equal to two explicit time steps on the fine grid. 
It may be noted also that the computed drag coefficient of the Korn airfoil at the shock-free design point is 
zero to four digits. The drag coefficient is also computed to be zero to four digits for subsonic flows over a 
variety of airfoils with lift coefficients in the range up to 1.0. Very little change is observed between solutions 
calculated on 80x16 and 160x32 meshes, providing a further confirmation of accuracy. 

The CUSP scheme produces very sharp shock waves in hypersonic flow, provided that care is taken to 
define the cell interface Mach number as the Mach number on the downwind side, so that downwind terms 
are perfectly canceled in supersonic flow. This is illustrated in Figures 10 and 11, which show the flow past a 
semicircular blunt body at Mach 8 and 20. It can be seen that quite rapid convergence, at a rate of the order 
of 0.9, continues to be obtained with the multigrid scheme in hypersonic flow. 


5.3 Steady flow calculations using characteristic splitting 

The remaining figures show the results of using characteristic splitting, with Roe linearization of the Jacobian 
matrices across the cell interfaces. The difference scheme was the symmetric limited positive (SLIP) scheme, 
applied to the differences of the characteristic variables. Let 


where A/ l+ ^ j = -A+i.j — Awi+$,j = tu,*+ij — Wij and is the Jacobian matrix calculated with Roe 

averaging, 

y/Pi+lj + y/PiJ 


for any quantity q. Then define 

= T ~ lAw i+i,j 
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where A i+ ^ j = TAT -1 , and apply the SLIP construction to each element of Ai> l+ i j separately. For each 
characteristic field the diffusion coefficient is taken as 

+ i = 2 A 

where A = |A| if |A| > e , + ^r) if |A| < e. The limiter defined by equation (18) was used, with q = 3, r = 

With this limiter the scheme is also precisely the JST scheme defined by equations (13) and (14). Therefore 
the scheme may be designated JSTR (Jameson-Schmidt-Turkel-Roe). 

Figures 12 and 13 show calculations for the same two cases as Figures 7 and 8, the RAE 2822 airfoil at Mach 

0.75 and 3° angle of attack, and the NACA 0012 airfoil at Mach 0.8 and 1.25° angle of attack, calculated on a 
very fine 320x64 meshes with O- topology. In each case the convergence history is shown for 200 cycles, while 
the pressure distribution is displayed after a sufficient number of cycles for its convergence. Only 25 cycles 
were needed for the RAE 2822 airfoil. Convergence was slower for the NACA 0012 airfoil, because additional 
cycles are needed to damp out a wave downstream of the very weak shock wave on the lower surface. Both 
calculations verify that the SLIP or equivalent JST schemes can resolve shocks with just one interior point if 
they are combined with characteristic splitting. 

Figures 14 and 15 show the performance of the scheme in hypersonic flow past a semicircular blunt body 
for the same conditions as Figures 10 and 11, Mach 8 and 20. These figures also exhibit the sharp discrete 
shocks which are obtained with the SLIP-JSTR construction. 500 multigrid cycles were used in each of these 
calculation. While the convergence is not as fast as the CUSP scheme, the multigrid scheme is still quite 
effective. 

Figures 16 and 17 show applications of the SLIP scheme with characteristic splitting to two airfoils which had 
previously been found to have non-unique solutions in calculations using the JST scheme with scalar diffusion 
[18]. The non-uniqueness is confirmed in these calculations, supporting the belief that these airfoils truly 
admit non-unique Euler solutions. The conditions under which the nonuniqueness was verified are identical 
with those found with the earlier scheme. In order to force the selection of one or other- of the solutions, it 
was sometimes necessary to start the calculation at a slightly higher or lower angle of attack, and then shift 
it by .05° after 200 cycles. This can be seen in the convergence histories. 

Finally, Figure 18 shows a three-dimensional Euler solution for the ONERA M6 wing at Mach 0.840 and 
an angle of attack of 3.06° calculated on a 192x32x48 mesh with 0-0 topology using the SLIP scheme with 
characteristic splitting. This again verifies the non-oscillatory character of the solution and sharp resolution 
of shock waves. 


6 Conclusion 

These numerical experiments confirm the theory of local extremum diminishing (LED) schemes, as it has been 
set forth in this paper. The following are the main conclusions of this study: 

1. The scalar diffusion that has been widely used can be significantly improved by the addition of a pressure 
term as defined in the CUSP formulation. Sharp discrete shocks are then obtained at high Mach numbers, 
and rapid multigrid convergence at all Mach numbers. 

2. The use of a split diffusive flux corresponding to the characteristic fields with Roe averaging improves the 
resolution of shocks in the transonic range, particularly when they are weak. 

3. The switched Jameson-Schmidt-Turkel (JST) scheme with the improved switch defined by Equations (13) 
and (14), and the equivalent symmetric limited positive (SLIP) scheme, defined by equations (15) and 
(18), are effective for steady state calculations in a wide Mach range. 

4. Corresponding symmetric limited positive (SLIP) and upstream limited positive (USLIP) schemes can be 
defined for both structured and unstructured meshes. 
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In the continuation of this paper the construction of numerical fluxes for the gas dynamic equations is 
further examined, and conditions are found under which steady discrete shock waves can contain a single 
interior point. It is shown that while characteristic decomposition is one way to achieve this property, it is not 
the only way, and equally sharp resolution can be obtained with less complex splittings. 
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Velocity 


6a: Pressure. 


6b: Density. 



6c: Velocity. 6d: Energy. 


Figure 6: Shock Tube Problem using SLIP Scheme with Pressure and Density Ratios of 10.0 and 8.0, 
respectively. Computed Results (+) are Compared with Theory ( — ) for 160 Equally Spaced Mesh Points 




/ 

i 


***>Q*I» 



7a: C p after 12 Cycles. 
Q = 1.1331, Cd = 0.0474. 



7b: C v after 50 Cycles. 
Ci — 1.1341, C d = 0.0475. 
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7c: Convergence. 7 d : Grid. 


Figure 7: RAE 2822. Mach 0.750, Angle of Attack 3°, 160x32 Mesh. 

CUSP Scheme. 


8a: C p after 12 Cycles. 
C, = 0.3672, C d = 0.0235. 


8b: C p after 50 Cycles. 
C, = 0.3688, C d = 0.0236 



8c: Convergence. 8d: Grid. 


Figure 8: NACA 0012. Mach 0.800, Angle of Attack 1.25°, 160x32 Mesh. 

CUSP Scheme. 


9a: C p after 12 Cycles. 
C, = 0.6317, C d = 0.0000. 


9b: C p after 50 Cycles. 
Ci = 0.6314, C d = 0.0000. 
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9c: Convergence. 



9d: Grid. 


Figure 9: KORN Airfoil. Mach 0.750, Angle of Attack 0°, 160x32 Mesh. 

CUSP Scheme. 



10a: C p on the Centerline in Front of the Body. 


10b: Convergence. 


Figure 10: Bluff Body. Mach 8, 160x64 Mesh 
CUSP Scheme. 



11a: C p on the Centerline in Front of the Body. lib: Convergence. 


Figure 11: Bluff Body. Mach 20, 160x64 Mesh. 
CUSP Scheme. 



12a: C v after 25 cycles. Woric 

C\ = 1.1203, Cd = 0.0457. 12b: Convergence. 

Figure 12: RAE 2822 with SLIP-JST Scheme with Characteristic Splitting. 
Mach 0.750, Angle of Attack 3°, 320x64 Mesh. 



13a: C p after 100 cycles. w,rt 

Ci = 0.3597, Cd — 0.0229. 13b: Convergence, 

Figure 13: NACA 0012 with SLIP-JST Scheme with Characteristic Splitting. 
Mach 0.800, Angle of Attack 1.25°, 320x64 Mesh. 
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14a: C p on the Centerline in Front of the Body. 


14b: Convergence. 


Figure 14: Bluff Body. Mach 8, 160x64 Mesh. 
SLIP-JST Scheme with Characteristic Splitting. 



15a: C p on the Centerline in Front of the Body. 



15b: Convergence 


Figure 15: Bluff Body. Mach 20, 160x64 Mesh. 
SLIP-JST Scheme with Characteristic Splitting. 
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10 a: L>p . 

Ci = 0.5797, C d = 0.0063. 


16b: Convergence. 




16d: Convergence. 


Figure 16: J-78 Airfoil: Non-Unique Solutions. 
SLIP Scheme with Characteristic Splitting. 
Mach 0.780, Angle of Attack -0.60°, 321x65 Mesh. 


i ' a: '-'p- 17b: Convergence. 

Ci = 0.5578, C d = 0.0050. 



17c: C p . 17d: Convergence. 

Ci = 0.6069, C d = 0.0010. 


Figure 17: GAW75-06-15 Airfoil: Non-Unique Solutions. 
SLIP Scheme with Characteristic Splitting. 

Mach 0.750, Angle of Attack -2.250°, 321x65 Mesh. 




18c: 50.00% Span. 18d: 68.75% Span. 

C, = 0.3261, C d = 0.0083. C, = 0.3190, C d = 0.0018. 


Figure 18: Onera M6 Wing. 

Mach 0.840, Angle of Attack 3.06°, 192x32x48 Mesh. 
C L = 0.3048, C D = 0.0125. 

SLIP Scheme with Characteristic Splitting. 
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